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A C++ INCARNATION OF ZERNIKE CIRCLE FUNCTIONS 


RICHARD J. MATHAR 


Abstract. An explicit CH— h library is provided which deals with Zernike 
Functions over the unit circle as the main subject. The implementation in¬ 
cludes basic means to evaluate the functions at points inside the unit circle 
and to convert the radial and azimuthal parameters to Noll’s index and vice 
versa. Advanced methods allow to expand products of Zernike Functions into 
sums of Zernike Functions, and to convert Zernike Functions to polynomials 
over the two Cartesian coordinates and vice versa. 


1. Basis Functions 


The Zernike circle functions Z are products of a radial polynomial R and an 
azimuthal sine or cosine function A\ 

(1) Z(™)(r,^)=R-(r)A„(^), 
where 

(2) r = 0 < r < 1 

is the distance to the origin of coordinates and 

(3) x = rcos(p] y = rsin(p 

define the azimuth angle ip. n is the degree of R, a non-negative integer, m is 
one of —n, —n + 2,..., n — 2, n, such that n — m is an even integer number. The 
normalization chosen here is [4] 


( 4 ) 


/ rR^{r)R'^,{r)dr = 5n,n'', / Am{‘p)A,n'{p)d(p = 

Jo Jo 


( 5 ) 

such that 

( 6 ) 


rdr / dipZ^(r,ip)Z^ (r,ip) = Sn,n'Sm,r. 


/2±M \ 

R”\r) ^ y2T+2(-l)("-l-l)/^(^„_^^|jrM2fi(- 


n — |ml n + ItoI , , n. 


and 

( 7 ) 


A (,n\ - / cos(m(^)/yi;;^, m > O; 

“ 1 sin(|m|(^)/y^, m<0. 
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where 


( 8 ) 


f 2, m = 0; 
1 , m ^ 0. 


Remark 1. This is sligthly different from the (more common) notation in my 
earlier representation [2, 3]; therefore some of the equations are reproduced here 
where factors need to be replaced. 

The product expansions of the azimuth functions are [2, §11.E]: 

( 9 ) 

^ ( ^y^\rn—rrl'\-^\7n—rn'\ -y/-^m+m ^: 777 ^ 0, 777- ^0, 

AmAm' = X I A_|^+|^/|| - Sgn(777 - |777'| )^_|^_ |^/11, TU > 0,777' < 0; 

Vemem'TT y j “ ^|m+m'|, 771 < 0, m' < 0. 

The product expansions of the radial functions are: 


( 10 ) 


n+n' 

K{r)R0:{r)= 


in.m.n' ,m' ,n" 


’K" ir), 


with projections [2, §11.E] 


1 3 


(11) 5ni,mi,n2,m2,n3,m3 / ^ (^)^' 

7o 


3 (r!,i-|mi|)/2 (ri2-|m2|)/2 (n3-|m3|)/2 

. ^IlK + i) E E E 

\ j=l SI—0 S2—0 S3—0 

xflH 


1 


2 + 77l + 772 + 773 — 2(si + S2 + S3) 


Sj I 


77 j 2Sj 


Sj J \(nj - lmjl)/2-SjJ' 


The class functions which are described in the next section answer the following 
questions: 

(1) Given the numerical coefficients Cn,m and Cn',m', what are the values of 


the linearization coefficients c„ 


in V r Y' r , ,7'-" 


(m') 




r^{m )c^ 

r," ^n" ,m" ^r)" 


(2) Given the numerical coefficients what are the values of the coefficients 

Cn,m in ~ '^n,Tn ^rL,mZn 

(3) Vice versa, given the numerical coefficients Cn,rm what are the values of the 

coefficients Cp,, in m = Sp q 


2. Implementation 

2.1. Points, Locations. Points in the unit circle and in the unit sphere are repre¬ 
sented by their 2 to 3 Gartesian coordinates and implemented in the classes Point2D 
and PointSD. The constructors expect the 2 or 3 Gartesian coordinates that fix the 
point in space. The Point2D object has a trivial method to extract the radial 
distance r and the azimuth angle ip of the circular coordinates from the x and y 
components. 
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2.2. Polynomials. A single term cx^ of a univariate polynomial is represented by 
an object of the class MonomiallD. A single term of a bivariate polynomial 

is represented by an object of the class Monomial2D. A single term cx^y'^z^ of a 
trivariate polynomial is represented by an object of the class MonomialSD. The 
constructors accept the exponents j, the exponents p and q and the exponents p, 
q and r respectively. The coefficient c is optional and set to unity if missing in the 
constructor. The useful operations within these classes are multiplication with or 
division through a constant value (which scales the coefficient), and multiplication 
with a single term of the same type (which essentially adds the exponents of the 
two factors). These operations are implemented by overloaded multiplication and 
division operators. There is a common method at () which evaluates the term given 
a value x on the line or a position specified by a Point2D or PointSD in the plane 
or in three-dimensional space. 

Univariate polynomials CjX^, bivariate polynomials ^ Cp^qX^y'^ and trivariate 
polynomials ^ are represented as vectors of the monomial objects in 

the classes PolynomiallD, Polynomial2D and Polynomial3D. Addition, subtraction 
and multiplications within each of them are closed operations and implemented by 
overloaded addition, subtraction and multiplication operators. Adding new terms 
is supported with overloaded += operators that accept a single monomial object or 
another polynomial object of the same dimension. 

There is a common method at () which evaluates these polynomials given a value 
X on the line or a position specified by a Point2D or PointSD. The function sums 
up the single-term components. 

A special variant of the univariate polynomials are the terminating Gaussian 
Hypergeometric Functions 2 F 1 (a, b; c; z) with integer parameter a < 0, two auxiliary 
parameters b and c, as a function of the variable t. These are implemented as a 
class Hypergeom21 derived from PolynomiallD. 


2.3. Zernike Radial Function. A radial function R^{r) is special case of the 
Hypergeometric Function constructed as in Equation (6) given the parameter n > 0 
and the parameter m (the latter being referred to only as |m|). It is represented by 
an object of the class ZernikeRadi, a subclass of Hypergeom21. A method in the 
class computes the g-factors of Equation (11) to support expansion of products of 
i?-functions in other i?-functions. 

The implementation is wider than actually needed here: an additional argument 
(which defaults to 2) supports use of the radial function in H > 2 dimensions, where 
D appears on the right hand sides of Eqs. (6) and (11) [3]. 


2.4. Zernike Azimuthal Function. A sine or cosine term of the form cAm is 
represented by an object of the class ZernikeCircAzi, which is constructed given 
the signed integer parameter m and an optional prefactor c—which is set to unity 
of missing. Scaling of the prefactor by a constant is implemented by overloaded 
multiplication and division. Evaluation at some explicit angle ip happens by calling 
the at() member function with an argument ip in units of radians. 

A collection (arithmetic sum) of the form is represented by an 

object of the class ZernikeCircAziVec, where each component is stored as an el¬ 
ement of the ZernikeCircAzi type. Adding new terms is achieved by using the 
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overloaded += operator. Scaling all terms (prefactors) at the same time is sup¬ 
ported by overloaded *= operation with a constant. An overloaded multiplication 
implements the arithmetic multiplication by means of Equation (9). 

2.5. Zernike Circle Function. A Zernike Circle Function cZ^'^ is represented 
by an object of the class ZernikeCirc which is initialized by the radial parameter 
n, the azimuthal parameter m and an optional prefactor c—which is set to unity 
if missing. Violation of the two constraints on m (parity and range) are silently 
caught by setting the prefactor to c = 0. Instead of the two parameters n and m, 
Noll’s index j >1 can also be used to define an object of ZernikeCirc. 

Scaling of the prefactor c is supported by overloaded *= and /= operators. 

Evaluation of the function at some point in the unit circle is done by calling 
the at() member function with an argument/location specified by a Point 2D ob¬ 
ject. The implementation constructs the factors i?™ and Am with objects of the 
ZernikeRadi and ZernikeCircAzi types and multiplies their values. A design 
decision is to set the values to zero outside the unit circle, r > !(!). 

An arithmetic sum of the form ^ Cn,mZ^'^ is represented by an object of the 
ZernikeCircVec class, which represents the components as members of a vector of 
ZernikeCirc’s. Adding terms to the sum is supported by overloaded += operators 
for additional ZernikeCirc or ZernikeCircVec terms. Scaling all the coefficients 
Cn,m is achieved by overloaded *= and /= operations with constant arguments. 

The most valuable functions of the implementation are: 

(1) The arithmetic product of two Zernike expansions is implemented by the 
overloaded operator *= which works with two factors of the ZernikeCirc or 
ZernikeCircVec type and produces a ZernikeCircVec object. This is an 
incarnation of [2, §11.E]; it starts with the product expansion of the A-terms 
provided the ZernikeCircAziVec class and collects the product expansion 
of the i?-terms by calling the g member functions in the ZernikeRadi class. 

(2) The conversion of a polynomial ^ cx^y'^ form into a Zernike expansion 
is supported by constructors in the ZernikeCircVec class that accept ar¬ 
guments of the Monomial2D or Polynomial2D type. This implements [2, 
§II.C]. 

(3) The conversion of a Zernike expansion ^ Cn,mZ^'^ into a polynomial ^ cx^y^ 
is supported by constructors in the Polynomial2D class that accept argu¬ 
ments of the ZernikeCirc or ZernikeCircVec type. This implements [2, 
§II.D]. 

2.6. Polynomial Fit. If the GNU Scientific Library is available [1], an additional 
class PowFit2D is introduced, which is fed with (constructed from) a list of Cartesian 
X and y coordinates and function values f{x,y) at these points. There are two 
simple formats of entering the data into the constructor, one that reads the x,y,f 
triples from an ASCII file, the other providing them as a list of PointSD points 
(interpreting the third coordinate as /). 

The major member function is the f it () function which constructs the ordinary 
least squares fit through the / values up to some total order of the fitting polynomial 
^ ax^y ’^. The limiting order p-P g is an argument of the f it (). This is implemented 
as another constructor for the Polynomial2D class that admits a PowFit2D object 
as its argument. Fitting of scattered data over the unit circle to Zernike Circle 
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Polynomials for some upper limit of the index n is then a matter of converting the 
fitting Polynomial2D to a ZernikeCircVec object in a final step. 

Note that this is a demonstration of library usage, but the strategy is inefficient. 
In practise, the Zernike coefficients would be obtained directly by solving the linear 
algebra in the (r, (p) coordinates. 
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Appendix A. Installation 

A.l. Compilation. The source code of roughly 4300 lines is available in the anc 
directory, licensed under the GNU General Public License. 

It is compiled with the GNU autotools [5] via 
autoreconf -i 

./configure —prefix=$HDME 
make 

make install 

This bundles the classes in a library libZernikeCirc. a and compiles the test 
routine tstZernikeCirc. It also searches for the gsl and gslcblas libraries and 
includes PowFit2D if these are found. 

If the autotools are not available, compilation in the conventional style of 
g++ -c -02 [A-Z]*.cxx 
Id -i -o libZernikeCirc.a *.o 

g++ -o tstZernikeCirc tstZernikeCirc.cxx -L. -1 ZernikeCirc 
is an alternative. Definition of the preprocessor variables HAVE_GSL_GSL_SF_H and 
HAVE_GSL_GSL_L1NALG_H and adding the flags -Igsl -Igslcblas must be done 
manually then, if applicable. 

If doxygen is available, the API documentation can be constructed in the html 
directory with 
make doc 

firefox html/index.html 

A.2. Numerical Tests. The test program can be run with 

tstZernikeCirc 

The test suite contains 

(1) a table of mappings of the two parameters (n, m) onto Noll’s index j to test 
the nollldx member function of ZernikeCirc; 

(2) a table of mappings of Noll’s index j to the two parameters (n, m) to test 
basically the inverse functionality in the constructor of ZernikeCirc; 

(3) a double loop over pairs of objects of ZernikeCirc to test that their prod¬ 
ucts ZernikeCircVec have the same value at some Point2D as expected 
from the product of the individual values; 
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(4) a loop over various terms of the cx^y’^ format to test that the conversion 
of a Monomial2D object into a ZernikeCircVec object keeps its value for 
some points Point2D scattered in the unit circle; 

(5) a loop over various terms of the cZ^'^ format to test that the conversion of 
a ZernikeCirc object into a Polynomial2D object keeps its value for some 
points Point2D scattered in the unit circle. 
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